#!/usr/bin/python -itt
import sys
import numpy as np

import logging
import logging.config

import os
scriptpath = os.path.dirname(os.path.realpath(__file__))
logging.config.fileConfig(scriptpath+'/logging.conf')
logger = logging.getLogger('absorbPhotons')


def _absorb_photons( thedict,
        dictname='photonData',
        absorption_key='wavelength',
        absorption_func_key = 'survival_prob_wavelength'):

    """ Delete photons from thedict[dictname]
    
    This function assumes to get a dict full of numpy arrays, photonData.
    The user should provide a function-like *absorption_func*, which returns a survival probability 
    for a photon according to its property *absorption_key*.
    Or short::
    
       survival_probability = absorption_func( photonData[absorption_key] )
    
    _absorb_photons() decides for each photon based on its survival_probability,
    whether it survives or not. After this step the absorbed photons are deleted 
    from all numpy arrays in *photonData*
    
    In a future version, just a boolean numpy array might be generated, which simply 
    stores, if a photon survived or not. Calculations, that follow after such an
    "unfinished" absorption have to take into account this fact. 
    I might also need numpy.ma.MaskedArray for this, I am not sure yet.
    """
    
    absorption_func = thedict['settings'][absorption_func_key]
    
    if absorption_func is None:
        logger.warning("_absorb_photons called without absorption_func. Nothing was absorbed.")
        return thedict
    
    logger.debug("_absorb_photons starting up with:")
    logger.debug("   var_name:%s %s" % (absorption_key,absorption_func_key))
    
    # get that particular variable of all photons, which are the reason for the absorption
    # e.g. key='wavelength'
    var = thedict[dictname][absorption_key]
    
    # calulate probability, for each photon to survive.
    probability = absorption_func(var)
    logger.debug("probability.shape %r" % probability.shape)
    
    # get as many random values between 0.0 and 1.0 as we have photons
    num_photons = thedict[dictname][absorption_key].shape[0]
    rndm = np.random.rand(num_photons)
    logger.debug("rndm.shape %r" % rndm.shape)
    
    # if the randomvalue is larger than the probability, the photon is absorbed
    # if the randomvalue is smaller, the photon survives
    survivors_ids = np.where(rndm < probability)[0]
    logger.debug("absorption key: %s" % absorption_key)
    logger.debug("num_photons=%d" % num_photons)
    logger.debug("photons not absorbed=%d --> %0.2f percent" % (len(survivors_ids), float(len(survivors_ids))/num_photons*100))

    # now I have to perform the actual absorption, i.e. 
    # I have to delete photons from the run
    # therefor I loop over run.d, which is a dict full of numpy arrays
    for key in thedict[dictname]:
        thedict[dictname][key] = thedict[dictname][key][survivors_ids]
    return thedict
    
def absorb_photons( run ):
    """ Apply two absorptions and one turn before reflector starts.
    
    This function is usually called after readcorsika.readCorsikaFile() like::
    
       run = readcorsika.readCorsikaFile( filename )
       run = absorb_photons( run )
       
    Withing this function photons are being "absorped" according to the wavelength dependend 
    absorptions within the mirrors, the cones, and in a sense by the limited PDE of the SiPMs.
    Then the main incoming direction of the photons is turned into zenith and then 
    some photons are deleted, which are too far away from the reflector to be reflected at all.
    """
    import turnIntoZenit
    
    thedict = absorb_photons._absorb_photons( thedict, 
        dictname='photonData',
        absorption_key='wavelength',
        absorption_func_key = 'survival_prob_wavelength')

    thedict = turnIntoZenit.turn_photons_into_zenith( thedict, 
            dictname='photonData', 
            pos_key='ground_dir', 
            dir_key='ground_pos',
            new_pos_key='ground_pos_turned', 
            new_dir_key='ground_dir_turned')

    thedict = absorb_photons._absorb_photons( thedict, 
            dictname='photonData',
            absorption_key='ground_pos_turned',
            absorption_func_key = 'is_in_inner_circle_checker')
    
    # spw --> survival probability vs wavelength
    spw = fact_si.interpolate_from_files( 
                    fact_si.Wavelength_dependency_description_files )
    data = run['photonData']
    data = _absorb_photons( data, absorption_key='wavelength', absorption_func=spw)
    data = turnIntoZenit.turn_photons_into_zenith( data )
    
    isInInnerCircle = IsInInnerCircle(210)
    data = _absorb_photons( data, absorption_key='ground_pos_turned', absorption_func=isInInnerCircle)
    
    return run
    
if __name__ == '__main__':

    from optparse import OptionParser
    
    parser = OptionParser()
    parser.add_option("-f", "--file", dest="filename",
        help="write report to FILE", metavar="FILE")
    #parser.add_option("-p", "--plot",
    #    action="store_true", dest="create_plots", default=False,
    #    help="create some informative debugging plots")
    options, args = parser.parse_args()
    if options.filename is None:
        print 'Please submit a filename. "', sys.argv[0],'-h" for help'
        sys.exit(1)

    import readcorsika
    run = readcorsika.readCorsikaFile( options.filename )
    
    run = absorb_photons(run)
